WL2_Growth_Climate
Brandie QC
2025-05-07
WL2 Growth + Climate
Add climate to best base model
Weekly climate as predictor of K
Summary of climate as predictor of Hmax and maybe delta
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(brms)## Loading required package: Rcpp
## Loading 'brms' package (version 2.21.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
##
## The following object is masked from 'package:stats':
##
## arlibrary(zoo) #rollmean##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericLoad the growth data
all_wl2_size_all_nobuffers <- read_csv("../output/WL2_Traits/WL2-2023_Size_Combined.csv")## Rows: 17336 Columns: 10
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (5): block, Genotype, pop.mf, parent.pop, survey.notes
## dbl (4): mf, rep, height.cm, long.leaf.cm
## date (1): survey_date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Pop Location Info
wl2_gowers_2023 <- read_csv("../output/Climate/Gowers_WL2.csv") %>%
select(parent.pop:GrwSsn_GD, Wtr_Year_GD) %>%
pivot_wider(names_from = TimePd, values_from = c(GrwSsn_GD, Wtr_Year_GD)) ## Rows: 46 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): parent.pop, elevation.group, TimePd
## dbl (9): elev_m, Lat, Long, GrwSsn_GD, GrwSsn_FLINT_GD, GrwSsn_BIOCLIM_GD, W...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(wl2_gowers_2023)## # A tibble: 6 × 9
## parent.pop elevation.group elev_m Lat Long GrwSsn_GD_Recent
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 YO4 High 2158. 37.8 -120. 0.216
## 2 WL2 High 2020. 38.8 -120. 0.242
## 3 BH Low 511. 37.4 -120. 0.264
## 4 SQ3 High 2373. 36.7 -119. 0.265
## 5 SQ2 Mid 1934. 36.7 -119. 0.270
## 6 SQ1 Mid 1921. 36.6 -119. 0.271
## # ℹ 3 more variables: GrwSsn_GD_Historical <dbl>, Wtr_Year_GD_Recent <dbl>,
## # Wtr_Year_GD_Historical <dbl>Merge
all_wl2_size_loc <- left_join(all_wl2_size_all_nobuffers, wl2_gowers_2023) %>%
mutate(day=if_else(survey_date=="2023-07-03" | survey_date=="2023-07-06", 0, #day 0 = pre-transplant
as.numeric(survey_date-min(survey_date)-7))) #days in the field starts on 7/11 (day all plants taken to the field) rather than the exact planting date ## Joining with `by = join_by(parent.pop)`head(all_wl2_size_loc)## # A tibble: 6 × 19
## survey_date block Genotype pop.mf parent.pop mf rep height.cm
## <date> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 2023-07-03 <NA> CP2_1_1 CP2_1 CP2 1 1 0.5
## 2 2023-07-03 <NA> CP2_1_2 CP2_1 CP2 1 2 0.7
## 3 2023-07-03 <NA> CP2_1_3 CP2_1 CP2 1 3 1.1
## 4 2023-07-03 <NA> CP2_1_4 CP2_1 CP2 1 4 0.8
## 5 2023-07-03 <NA> CP2_1_5 CP2_1 CP2 1 5 0.9
## 6 2023-07-03 <NA> CP2_1_6 CP2_1 CP2 1 6 1
## # ℹ 11 more variables: long.leaf.cm <dbl>, survey.notes <chr>,
## # elevation.group <chr>, elev_m <dbl>, Lat <dbl>, Long <dbl>,
## # GrwSsn_GD_Recent <dbl>, GrwSsn_GD_Historical <dbl>,
## # Wtr_Year_GD_Recent <dbl>, Wtr_Year_GD_Historical <dbl>, day <dbl>Data set for the model (remove WV and shrinkage months)
all_wl2_size_for_models <- all_wl2_size_loc %>%
filter(parent.pop!="WV") %>%
group_by(Genotype) %>%
mutate(height_next = lead(height.cm, order_by = survey_date),
height_diff=height_next-height.cm,
shrink_day=if_else(height_diff>-2, NA, day),
min_shrink_day=min(shrink_day, na.rm = TRUE)) %>% #deals with problem of if something shrinks and then starts to grow again
filter(day<min_shrink_day | min_shrink_day==Inf) %>%
ungroup() %>%
drop_na(height.cm)## Warning: There were 1392 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `min_shrink_day = min(shrink_day, na.rm = TRUE)`.
## ℹ In group 1: `Genotype = "BH_1_1"`.
## Caused by warning in `min()`:
## ! no non-missing arguments to min; returning Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1391 remaining warnings.Read in soil temp
temp_WL2 <- read_csv("../input/WL2_Data/WL2_2022_2023_iButton_Data_Corrected.csv") %>%
select(-`...3`) %>%
mutate(Date_Time = mdy_hm(Date_Time)) %>%
filter(Date_Time > ymd("2023-07-06"))## New names:
## Rows: 14253 Columns: 4
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): Bed, Date_Time, ...3 dbl (1): SoilTemp
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...3`head(temp_WL2)## # A tibble: 6 × 3
## Bed Date_Time SoilTemp
## <chr> <dttm> <dbl>
## 1 A_1 2023-07-13 14:46:00 27.5
## 2 A_1 2023-07-13 15:46:00 32
## 3 A_1 2023-07-13 16:46:00 31.5
## 4 A_1 2023-07-13 17:46:00 28.5
## 5 A_1 2023-07-13 18:46:00 25
## 6 A_1 2023-07-13 19:46:00 22skimr::skim(temp_WL2)| Name | temp_WL2 |
| Number of rows | 10240 |
| Number of columns | 3 |
| _______________________ | |
| Column type frequency: | |
| character | 1 |
| numeric | 1 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| Bed | 0 | 1 | 3 | 3 | 0 | 5 | 0 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| SoilTemp | 0 | 1 | 16.22 | 8.51 | 1.5 | 10.5 | 14 | 20 | 50.5 | ▅▇▂▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| Date_Time | 0 | 1 | 2023-07-13 14:40:00 | 2023-10-06 21:48:00 | 2023-08-25 06:14:00 | 10240 |
Average across beds.
temp_WL2_summary <- temp_WL2 %>%
group_by(Date_Time) %>%
summarise(AvgSoilTemp=mean(SoilTemp, na.rm=TRUE)) %>% #avg across beds
mutate(Date=as.Date(Date_Time)) %>%
filter(Date != min(Date), Date != max(Date)) %>% # trim ragged ends
group_by(Date) %>% #summmarise hourly data by date
summarize(
min_temp_d = min(AvgSoilTemp),
max_temp_d = max(AvgSoilTemp),
mean_temp_d = mean(AvgSoilTemp)
) %>%
mutate(
across(ends_with("temp_d"), \(x) rollmean(x, k = 7, align = "right", fill = NA), .names="{.col}1_7"),
across(ends_with("temp_d"), \(x) rollmean(x, k = 13, align = "right", fill = NA), .names="{.col}1_13") # 13 so I can get the first survey date
)
temp_WL2_summary## # A tibble: 84 × 10
## Date min_temp_d max_temp_d mean_temp_d min_temp_d1_7 max_temp_d1_7
## <date> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2023-07-14 8.5 31 17.0 NA NA
## 2 2023-07-15 11 35 18.7 NA NA
## 3 2023-07-16 11 36 19.1 NA NA
## 4 2023-07-17 14 36 20.5 NA NA
## 5 2023-07-18 13 36 19.8 NA NA
## 6 2023-07-19 12.5 36 19.3 NA NA
## 7 2023-07-20 11.5 45 21.4 11.6 36.4
## 8 2023-07-21 11 47.5 21.9 12 38.8
## 9 2023-07-22 11.5 48.5 22.5 12.1 40.7
## 10 2023-07-23 14 48.5 23.4 12.5 42.5
## # ℹ 74 more rows
## # ℹ 4 more variables: mean_temp_d1_7 <dbl>, min_temp_d1_13 <dbl>,
## # max_temp_d1_13 <dbl>, mean_temp_d1_13 <dbl>Read in soil moisture
moisture_WL2 <- read_csv("../input/WL2_Data/WL2_2023_Bed_C_Soil_Moisture_Corrected.csv") %>%
mutate(Date_Time = mdy_hm(Date_Time)) ## Rows: 2209 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): Date_Time
## dbl (5): Port_1, Port_2, Port_3, Port_4, Port_5
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(moisture_WL2)## # A tibble: 6 × 6
## Date_Time Port_1 Port_2 Port_3 Port_4 Port_5
## <dttm> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2023-07-20 12:00:00 0.071 0.141 0.075 0.055 0.092
## 2 2023-07-20 13:00:00 0.07 0.142 0.077 0.056 0.094
## 3 2023-07-20 14:00:00 0.074 0.142 0.079 0.056 0.094
## 4 2023-07-20 15:00:00 0.076 0.141 0.081 0.055 0.094
## 5 2023-07-20 16:00:00 0.077 0.14 0.082 0.054 0.093
## 6 2023-07-20 17:00:00 0.077 0.136 0.081 0.052 0.091skimr::skim(moisture_WL2)| Name | moisture_WL2 |
| Number of rows | 2209 |
| Number of columns | 6 |
| _______________________ | |
| Column type frequency: | |
| numeric | 5 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| Port_1 | 0 | 1 | 0.13 | 0.06 | 0.00 | 0.08 | 0.14 | 0.17 | 0.36 | ▆▇▇▁▁ |
| Port_2 | 0 | 1 | 0.10 | 0.02 | 0.04 | 0.08 | 0.10 | 0.12 | 0.16 | ▁▅▇▇▂ |
| Port_3 | 0 | 1 | 0.10 | 0.05 | -0.01 | 0.06 | 0.12 | 0.13 | 0.21 | ▃▂▅▇▁ |
| Port_4 | 0 | 1 | 0.03 | 0.02 | -0.05 | 0.01 | 0.03 | 0.05 | 0.10 | ▁▂▇▆▁ |
| Port_5 | 0 | 1 | 0.08 | 0.04 | -0.03 | 0.06 | 0.09 | 0.10 | 0.19 | ▂▂▇▅▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| Date_Time | 0 | 1 | 2023-07-20 12:00:00 | 2023-10-20 12:00:00 | 2023-09-04 12:00:00 | 2209 |
Will average across the ports
moisture_WL2_summary <- moisture_WL2 %>%
mutate(Port_1=if_else(Port_1<0, 0, Port_1),
Port_2=if_else(Port_2<0, 0, Port_2),
Port_3=if_else(Port_3<0, 0, Port_3),
Port_4=if_else(Port_4<0, 0, Port_4),
Port_5=if_else(Port_5<0, 0, Port_5)) %>% #convert negatives to zeros (negatives result from extremely dry soil --> air pockets)
rowwise() %>%
mutate(s_moisture = mean(c_across(-Date_Time)) ) %>%
select(Date_Time, s_moisture) %>%
mutate(Date=as.Date(Date_Time)) %>%
group_by(Date) %>%
summarize(
mean_moisture_d = mean(s_moisture)
) %>%
mutate(
s_moisture_1_7 = rollmean(mean_moisture_d, k = 7, align = "right", fill = "extend"),
s_moisture_1_13 = rollmean(mean_moisture_d, k = 13, align = "right", fill = "extend")
)
moisture_WL2_summary## # A tibble: 93 × 4
## Date mean_moisture_d s_moisture_1_7 s_moisture_1_13
## <date> <dbl> <dbl> <dbl>
## 1 2023-07-20 0.0853 0.0714 0.0571
## 2 2023-07-21 0.0809 0.0714 0.0571
## 3 2023-07-22 0.0754 0.0714 0.0571
## 4 2023-07-23 0.0711 0.0714 0.0571
## 5 2023-07-24 0.0668 0.0714 0.0571
## 6 2023-07-25 0.0622 0.0714 0.0571
## 7 2023-07-26 0.0581 0.0714 0.0571
## 8 2023-07-27 0.0537 0.0669 0.0571
## 9 2023-07-28 0.0484 0.0622 0.0571
## 10 2023-07-29 0.0436 0.0577 0.0571
## # ℹ 83 more rowsadd climate data
all_wl2_size_for_models_climate <- all_wl2_size_for_models %>%
left_join(temp_WL2_summary, by = c("survey_date" = "Date")) %>%
left_join(moisture_WL2_summary, by = c("survey_date" = "Date")) %>%
filter(day>0) #remove day pre-field size b/c no climate data
all_wl2_size_for_models_climate %>%
ggplot(aes(group=Genotype, x=survey_date, y=height.cm, col=elev_m)) +
scale_colour_gradient(low = "#F5A540", high = "#0043F0") +
geom_line() + facet_wrap(~parent.pop, scales = "free")
#without day 0, min is different...Priors
prior1_wl2 <- c(set_prior("normal(5,3)", nlpar="Hmin"),
set_prior("normal(60,15)", nlpar="Hmax"),
set_prior("gamma(2,2)", nlpar = "k", lb = 0), #gamma constrains to be > 0
set_prior("gamma(20,3)", nlpar = "delta", lb = 0)
)
#first number = avg, second number = the space around that average to search Base Models
#best base model:
base_f1 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
Hmin + k + delta ~ 1,
Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)
base_f2 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
k + delta ~ 1,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)Base Model Fits
Base model 1 - Hmax vary by pop
base_fit1 <- brm(formula=base_f1, data=all_wl2_size_for_models_climate, prior=prior1_wl2, iter = 4000, cores=4,
control = list(max_treedepth = 12))## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsaveRDS(base_fit1, file = "brm_object_wl2_base1_rptedmeas.rds")
summary(base_fit1)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## Hmin ~ 1
## k ~ 1
## delta ~ 1
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_climate (Number of observations: 6038)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmax_Intercept) 7.33 1.49 5.39 11.22 1.00 3514 2939
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmax_Intercept) 11.47 2.97 7.38 18.85 1.00 2458 2860
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Hmin_Intercept 2.64 0.07 2.50 2.77 1.00 8880 6398
## k_Intercept 0.41 0.17 0.13 0.77 1.00 3610 2851
## delta_Intercept 0.59 0.03 0.54 0.64 1.00 3702 3209
## Hmax_Intercept 9.32 2.93 4.49 15.99 1.00 1093 2238
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.07 0.01 1.05 1.09 1.00 10226 6189
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(base_fit1, nvariables = 3, ask=FALSE)
pred.df <- expand_grid(day=min(all_wl2_size_for_models_climate$day):max(all_wl2_size_for_models_climate$day),
Genotype=unique(all_wl2_size_for_models_climate$Genotype)) %>%
separate(Genotype, c("parent.pop", "mf","rep"), remove = FALSE)
pred.df2 <- all_wl2_size_for_models_climate %>% select(day, parent.pop, Genotype)
base_fit1.predictions <- cbind(pred.df2, prediction=predict(base_fit1, newdata = pred.df2)[,"Estimate"]) %>%
full_join(all_wl2_size_for_models_climate, by=c("day", "parent.pop", "Genotype"))
base_fit1.predictions %>% ggplot(aes(x=day)) +
geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
geom_line(aes(group=Genotype, y=prediction, color=parent.pop)) +
facet_wrap(~parent.pop, scales = "free")
ggsave("../output/WL2_GrowthClim_Predictions_Base1_Hmax_RptedMeas.png", width = 14, height = 8, units = "in")Base model 2 - Hmin + Hmax vary by pop
base_fit2 <- brm(formula=base_f2, data=all_wl2_size_for_models_climate, prior=prior1_wl2, iter = 4000, cores=4,
control = list(max_treedepth = 12))## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsaveRDS(base_fit2, file = "brm_object_wl2_base2_rptedmeas.rds")
summary(base_fit2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## k ~ 1
## delta ~ 1
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_climate (Number of observations: 6038)
## Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
## total post-warmup draws = 8000
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.86 0.05 1.76 1.96 1.00 5146 6543
## sd(Hmax_Intercept) 5.19 0.41 4.49 6.08 1.00 2519 3973
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.39 0.25 1.00 1.95 1.00 3936 5209
## sd(Hmax_Intercept) 8.40 1.43 6.12 11.67 1.00 4214 4876
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## k_Intercept 0.83 0.10 0.64 1.02 1.00 2489 4084
## delta_Intercept 1.03 0.05 0.95 1.13 1.00 2279 3882
## Hmin_Intercept 3.24 0.31 2.64 3.87 1.00 2850 4409
## Hmax_Intercept 7.01 1.85 3.48 10.83 1.00 1980 3037
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.84 0.01 0.82 0.86 1.00 7479 6494
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(base_fit2, nvariables = 3, ask=FALSE)
pred.df <- expand_grid(day=min(all_wl2_size_for_models_climate$day):max(all_wl2_size_for_models_climate$day),
Genotype=unique(all_wl2_size_for_models_climate$Genotype)) %>%
separate(Genotype, c("parent.pop", "mf","rep"), remove = FALSE)
pred.df2 <- all_wl2_size_for_models_climate %>% select(day, parent.pop, Genotype)
base_fit2.predictions <- cbind(pred.df2, prediction=predict(base_fit2, newdata = pred.df2)[,"Estimate"]) %>%
full_join(all_wl2_size_for_models_climate, by=c("day", "parent.pop", "Genotype"))
base_fit2.predictions %>% ggplot(aes(x=day)) +
geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
geom_line(aes(group=Genotype, y=prediction, color=parent.pop)) +
facet_wrap(~parent.pop, scales = "free")
ggsave("../output/WL2_GrowthClim_Predictions_Base2_HmaxHmin_RptedMeas.png", width = 14, height = 8, units = "in")Model comparison
base_fit1 <- add_criterion(base_fit1, "loo")## Warning: Found 93 observations with a pareto_k > 0.7 in model 'base_fit1'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.base_fit2 <- add_criterion(base_fit2, "loo")## Warning: Found 308 observations with a pareto_k > 0.7 in model 'base_fit2'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.loo_compare(base_fit1, base_fit2)## elpd_diff se_diff
## base_fit2 0.0 0.0
## base_fit1 -1000.5 100.9#model with Hmin is better Temp Models
temp_f2 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ min_temp_d1_7,
Hmin + Hmax ~ (1|parent.pop)+ (1|Genotype),
nl=TRUE)
temp_f3 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ mean_temp_d1_7,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)
temp_f4 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ max_temp_d1_7,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)
temp_f5 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ min_temp_d1_13,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)
temp_f6 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ mean_temp_d1_13,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)
temp_f7 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ max_temp_d1_13,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)Temp Model Fits
Min D1-7 (fit2)
all_wl2_size_for_models_temp_fit2 <- all_wl2_size_for_models_climate %>% drop_na(min_temp_d1_7)
temp_fit2 <- brm(formula=temp_f2, data=all_wl2_size_for_models_temp_fit2, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start sampling#Warning: There were 257 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: Examine the pairs() plot to diagnose sampling problems
summary(temp_fit2) #some ESS kind of low ## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ min_temp_d1_7
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_temp_fit2 (Number of observations: 5403)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.78 0.05 1.68 1.88 1.00 4557 6223
## sd(Hmax_Intercept) 6.98 0.95 5.50 9.14 1.00 2483 3873
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.48 0.26 1.07 2.07 1.00 3818 5009
## sd(Hmax_Intercept) 10.60 2.12 7.28 15.54 1.00 3015 4361
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.08 0.05 0.99 1.17 1.00 2516 5128
## k_Intercept 0.62 0.10 0.43 0.82 1.00 2347 3975
## k_min_temp_d1_7 0.00 0.00 0.00 0.01 1.00 8888 5563
## Hmin_Intercept 3.28 0.32 2.65 3.91 1.00 2207 3628
## Hmax_Intercept 8.41 2.49 3.87 13.69 1.00 1339 2854
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.76 0.01 0.74 0.77 1.00 7138 6042
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(temp_fit2, nvariables = 3, ask=FALSE) #some Hmax chain weirdness 
Mean D1-7 (fit3)
all_wl2_size_for_models_temp_fit3 <- all_wl2_size_for_models_climate %>% drop_na(mean_temp_d1_7)
temp_fit3 <- brm(formula=temp_f3, data=all_wl2_size_for_models_temp_fit3, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(temp_fit3)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ mean_temp_d1_7
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_temp_fit3 (Number of observations: 5403)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.77 0.05 1.68 1.87 1.00 4391 5105
## sd(Hmax_Intercept) 6.89 0.92 5.46 9.07 1.00 2181 3419
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.47 0.25 1.06 2.05 1.00 3683 5031
## sd(Hmax_Intercept) 10.44 2.08 7.26 15.27 1.00 3310 4698
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.13 0.05 1.04 1.22 1.00 3105 5145
## k_Intercept 0.59 0.09 0.41 0.77 1.00 2143 3568
## k_mean_temp_d1_7 0.01 0.00 0.00 0.01 1.00 5490 4379
## Hmin_Intercept 3.26 0.32 2.64 3.89 1.00 2274 3285
## Hmax_Intercept 8.36 2.43 3.93 13.54 1.00 1521 2778
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.75 0.01 0.74 0.77 1.00 5932 6072
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(temp_fit3, nvariables = 3, ask=FALSE)
#temp_fit3.predictions <- cbind(pred.df, prediction=predict(temp_fit3, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit3.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)Max D1-7 (fit4)
all_wl2_size_for_models_temp_fit4 <- all_wl2_size_for_models_climate %>% drop_na(max_temp_d1_7)
temp_fit4 <- brm(formula=temp_f4, data=all_wl2_size_for_models_temp_fit4, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(temp_fit4)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ max_temp_d1_7
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_temp_fit4 (Number of observations: 5403)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.80 0.05 1.71 1.90 1.00 4281 5523
## sd(Hmax_Intercept) 5.63 0.64 4.62 7.12 1.00 2106 3038
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.37 0.24 0.99 1.91 1.00 4333 5168
## sd(Hmax_Intercept) 8.74 1.65 6.14 12.60 1.00 3431 4947
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.06 0.04 0.98 1.14 1.00 2880 4393
## k_Intercept 0.74 0.11 0.53 0.95 1.00 2008 3185
## k_max_temp_d1_7 0.00 0.00 0.00 0.00 1.00 2687 4353
## Hmin_Intercept 3.20 0.31 2.59 3.80 1.00 2968 4317
## Hmax_Intercept 7.32 1.99 3.67 11.46 1.00 1788 3177
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.75 0.01 0.74 0.77 1.00 5457 5091
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(temp_fit4, nvariables = 3, ask=FALSE)
#temp_fit4.predictions <- cbind(pred.df, prediction=predict(temp_fit4, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit4.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)Min D1-13 (fit5)
all_wl2_size_for_models_temp_fit5 <- all_wl2_size_for_models_climate %>% drop_na(min_temp_d1_13)
temp_fit5 <- brm(formula=temp_f5, data=all_wl2_size_for_models_temp_fit5, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(temp_fit5)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ min_temp_d1_13
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_temp_fit5 (Number of observations: 5403)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.77 0.05 1.68 1.87 1.00 3716 5460
## sd(Hmax_Intercept) 7.40 1.12 5.69 9.98 1.00 2465 3502
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.51 0.26 1.08 2.08 1.00 3869 4435
## sd(Hmax_Intercept) 11.12 2.32 7.53 16.48 1.00 3028 3924
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.12 0.05 1.02 1.23 1.00 2682 4907
## k_Intercept 0.57 0.10 0.39 0.77 1.00 2386 3483
## k_min_temp_d1_13 0.01 0.00 0.00 0.01 1.00 6306 4789
## Hmin_Intercept 3.28 0.33 2.63 3.95 1.00 2372 3972
## Hmax_Intercept 8.92 2.57 4.25 14.59 1.00 1517 2538
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.76 0.01 0.74 0.77 1.00 5861 5796
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(temp_fit5, nvariables = 3, ask=FALSE)
#temp_fit5.predictions <- cbind(pred.df, prediction=predict(temp_fit5, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit5.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)Mean D1-13 (fit6)
all_wl2_size_for_models_temp_fit6 <- all_wl2_size_for_models_climate %>% drop_na(mean_temp_d1_13)
temp_fit6 <- brm(formula=temp_f6, data=all_wl2_size_for_models_temp_fit6, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(temp_fit6)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ mean_temp_d1_13
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_temp_fit6 (Number of observations: 5403)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.78 0.05 1.68 1.87 1.00 4587 5397
## sd(Hmax_Intercept) 7.02 0.97 5.51 9.29 1.00 2045 3978
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.49 0.26 1.08 2.07 1.00 3420 4655
## sd(Hmax_Intercept) 10.60 2.14 7.26 15.54 1.00 2800 4449
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.11 0.05 1.02 1.22 1.00 3031 4472
## k_Intercept 0.59 0.10 0.41 0.79 1.00 1907 3753
## k_mean_temp_d1_13 0.00 0.00 0.00 0.01 1.00 7555 4460
## Hmin_Intercept 3.27 0.33 2.64 3.94 1.00 2278 3673
## Hmax_Intercept 8.51 2.48 4.04 13.78 1.00 1544 2947
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.76 0.01 0.74 0.77 1.00 6782 5845
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(temp_fit6, nvariables = 3, ask=FALSE)
#temp_fit6.predictions <- cbind(pred.df, prediction=predict(temp_fit6, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit6.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)Max D1-13 (fit7)
all_wl2_size_for_models_temp_fit7 <- all_wl2_size_for_models_climate %>% drop_na(max_temp_d1_13)
temp_fit7 <- brm(formula=temp_f7, data=all_wl2_size_for_models_temp_fit7, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(temp_fit7)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ max_temp_d1_13
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_temp_fit7 (Number of observations: 5403)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.79 0.05 1.69 1.88 1.00 4632 5990
## sd(Hmax_Intercept) 6.45 0.80 5.19 8.28 1.00 2471 3549
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.44 0.26 1.03 2.03 1.00 3958 4880
## sd(Hmax_Intercept) 9.87 1.95 6.83 14.47 1.00 3535 5059
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.06 0.04 0.98 1.14 1.00 3159 4830
## k_Intercept 0.66 0.10 0.47 0.86 1.00 2341 3512
## k_max_temp_d1_13 0.00 0.00 0.00 0.00 1.00 5504 4237
## Hmin_Intercept 3.25 0.32 2.64 3.89 1.00 2790 3996
## Hmax_Intercept 7.91 2.30 3.54 12.76 1.00 1490 2808
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.76 0.01 0.74 0.77 1.00 6152 6133
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(temp_fit7, nvariables = 3, ask=FALSE)
#temp_fit7.predictions <- cbind(pred.df, prediction=predict(temp_fit7, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#temp_fit7.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)save.image("../output/BQC_WL2_Growth_Climate.Rdata")Moisture Models
moisture_f2 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ s_moisture_1_7,
Hmin + Hmax ~ (1|parent.pop) + (1|Genotype),
nl=TRUE)
moisture_f3 <- brmsformula(height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta),
delta ~ 1,
k ~ s_moisture_1_13,
Hmin + Hmax ~ (1|parent.pop) ,
nl=TRUE)Moisture Model Fits
Min D1-7 (fit2)
all_wl2_size_for_models_moisture_fit2 <- all_wl2_size_for_models_climate %>% drop_na(s_moisture_1_7)
moisture_fit2 <- brm(formula=moisture_f2, data=all_wl2_size_for_models_moisture_fit2, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(moisture_fit2)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ s_moisture_1_7
## Hmin ~ (1 | parent.pop) + (1 | Genotype)
## Hmax ~ (1 | parent.pop) + (1 | Genotype)
## Data: all_wl2_size_for_models_moisture_fit2 (Number of observations: 6038)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~Genotype (Number of levels: 1035)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.86 0.05 1.76 1.96 1.00 4761 4988
## sd(Hmax_Intercept) 5.20 0.41 4.51 6.14 1.00 2774 3997
##
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 1.39 0.25 0.99 1.94 1.00 4104 4462
## sd(Hmax_Intercept) 8.44 1.43 6.08 11.71 1.00 4223 5303
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 1.03 0.04 0.95 1.12 1.00 2257 4565
## k_Intercept 0.81 0.09 0.62 0.99 1.00 2609 3720
## k_s_moisture_1_7 0.13 0.08 0.02 0.31 1.00 8524 4542
## Hmin_Intercept 3.25 0.31 2.64 3.89 1.00 2939 4300
## Hmax_Intercept 7.08 1.89 3.48 10.95 1.00 1889 3346
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.84 0.01 0.82 0.86 1.00 5407 5643
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(moisture_fit2, nvariables = 3, ask=FALSE)
#moisture_fit2.predictions <- cbind(pred.df, prediction=predict(moisture_fit2, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))
#
#moisture_fit2.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)Mean D1-13 (fit3)
all_wl2_size_for_models_moisture_fit3 <- all_wl2_size_for_models_climate %>% drop_na(s_moisture_1_13)
moisture_fit3 <- brm(formula=moisture_f3, data=all_wl2_size_for_models_moisture_fit3, prior=prior1_wl2, iter = 3800, cores=4)## Compiling Stan program...## Trying to compile a simple C file## Running /usr/lib/R/bin/R CMD SHLIB foo.c
## using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
## gcc -I"/usr/share/R/include" -DNDEBUG -I"/usr/lib/R/site-library/Rcpp/include/" -I"/usr/lib/R/site-library/RcppEigen/include/" -I"/usr/lib/R/site-library/RcppEigen/include/unsupported" -I"/usr/lib/R/site-library/BH/include" -I"/usr/lib/R/site-library/StanHeaders/include/src/" -I"/usr/lib/R/site-library/StanHeaders/include/" -I"/usr/lib/R/site-library/RcppParallel/include/" -I"/usr/lib/R/site-library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -fpic -g -O2 -ffile-prefix-map=/build/r-base-TYDrW1/r-base-4.5.0=. -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -c foo.c -o foo.o
## In file included from /usr/lib/R/site-library/RcppEigen/include/Eigen/Core:19,
## from /usr/lib/R/site-library/RcppEigen/include/Eigen/Dense:1,
## from /usr/lib/R/site-library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## /usr/lib/R/site-library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [/usr/lib/R/etc/Makeconf:202: foo.o] Error 1## Start samplingsummary(moisture_fit3)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: height.cm ~ Hmax - (Hmax - Hmin) * exp(-(k/100 * day)^delta)
## delta ~ 1
## k ~ s_moisture_1_13
## Hmin ~ (1 | parent.pop)
## Hmax ~ (1 | parent.pop)
## Data: all_wl2_size_for_models_moisture_fit3 (Number of observations: 6038)
## Draws: 4 chains, each with iter = 3800; warmup = 1900; thin = 1;
## total post-warmup draws = 7600
##
## Multilevel Hyperparameters:
## ~parent.pop (Number of levels: 22)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Hmin_Intercept) 0.97 0.24 0.59 1.52 1.00 1946 2974
## sd(Hmax_Intercept) 5.79 1.01 4.22 8.06 1.00 2056 2872
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## delta_Intercept 0.95 0.13 0.74 1.23 1.00 2043 2328
## k_Intercept 1.93 0.28 1.31 2.41 1.00 5820 3273
## k_s_moisture_1_13 0.70 0.45 0.10 1.80 1.00 9792 4541
## Hmin_Intercept 2.89 0.28 2.37 3.47 1.00 2176 2590
## Hmax_Intercept 5.72 1.26 3.21 8.35 1.01 732 1282
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 2.49 0.02 2.44 2.53 1.00 13345 5406
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(moisture_fit3, nvariables = 3, ask=FALSE)
#moisture_fit3.predictions <- cbind(pred.df, prediction=predict(moisture_fit3, newdata = pred.df)[,"Estimate"]) %>%
# full_join(all_wl2_size_for_models, by=c("day", "parent.pop"))#
#moisture_fit3.predictions %>% ggplot(aes(x=day)) +
# geom_line(aes(group=Genotype, y=height.cm, color=parent.pop), alpha=.2) +
# geom_line(aes(y=prediction, color=parent.pop)) +
# facet_wrap(~parent.pop)Temp Model Comparison
#base_fit1 <- add_criterion(base_fit1, "loo")
#base_fit2 <- add_criterion(base_fit2, "loo")
temp_fit2 <- add_criterion(temp_fit2, "loo")## Warning: Found 331 observations with a pareto_k > 0.7 in model 'temp_fit2'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.temp_fit3 <- add_criterion(temp_fit3, "loo")## Warning: Found 319 observations with a pareto_k > 0.7 in model 'temp_fit3'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.temp_fit4 <- add_criterion(temp_fit4, "loo")## Warning: Found 324 observations with a pareto_k > 0.7 in model 'temp_fit4'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.temp_fit5 <- add_criterion(temp_fit5, "loo") ## Warning: Found 328 observations with a pareto_k > 0.7 in model 'temp_fit5'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.temp_fit6 <- add_criterion(temp_fit6, "loo") ## Warning: Found 340 observations with a pareto_k > 0.7 in model 'temp_fit6'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.temp_fit7 <- add_criterion(temp_fit7, "loo") ## Warning: Found 324 observations with a pareto_k > 0.7 in model 'temp_fit7'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.loo_compare(temp_fit2, temp_fit3, temp_fit4, temp_fit5, temp_fit6, temp_fit7)## elpd_diff se_diff
## temp_fit4 0.0 0.0
## temp_fit3 -1.7 7.2
## temp_fit7 -5.2 10.1
## temp_fit6 -5.6 10.3
## temp_fit5 -6.7 13.3
## temp_fit2 -13.0 12.3Moisture Model Comparison
#base_fit1 <- add_criterion(base_fit1, "loo")
#base_fit2 <- add_criterion(base_fit2, "loo")
moisture_fit2 <- add_criterion(moisture_fit2, "loo")## Warning: Found 291 observations with a pareto_k > 0.7 in model 'moisture_fit2'.
## We recommend to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.moisture_fit3 <- add_criterion(moisture_fit3, "loo")
loo_compare(moisture_fit2, moisture_fit3)## elpd_diff se_diff
## moisture_fit2 0.0 0.0
## moisture_fit3 -5670.1 179.1save.image("../output/BQC_WL2_Growth_Climate.Rdata")